Models of Strong Interaction in Flat-Band Graphene Nanoribbons: Magnetic 

Quantum Crystals 
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Graphene based nanostructures exhibit flat electronic energy bands in their single-particle spec- 
trum. We consider interacting electrons in flat bands of zig-zag nanoribbons. We present a protocol 
for flat-band projection that yields interaction-only tight-binding models. We argue that, at low 
densities, flat bands can delocalize single-particle basis states to support ferromagnetic quantum 
crystal ground states. 
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I. INTRODUCTION 

Graphene based structures offer unique opportuni- 
ties to engineer electronic band structure by shape 
alonei^ Infinite graphene sheets exhibit a conic spec- 
trum but finite sized graphene nanostructures yield a 
surprisingly broad array of interesting band features. 
A subset of graphene nanostructures reveal flat bands. 
Theoretical work shows that flat bands can be found, 
e.g., at the edges of two-dimensional graphene^ in 
one-dimensional graphene nanoribbons,— ~— hydrogenated 
graphene nanoribbons,— graphene dots,— and graphene 
antidots^ 

Electrons in flat kinetic energy bands pose challeng- 
ing theoretical problems. The absence of any dispersion 
leaves the Coulomb interaction to govern the low energy 
physics. Many common approximations fail in the ex- 
treme flat-band limit. A single flat band cannot lead to 
intra-band screening as in ordinary Fermi liquids, e.g., 
two-dimensional graphene sheets.— Magnetic properties 
in bulk graphene in particular occur in a regime where 
large screening effects (allowed by a dispersive kinetic 
energy) minimize the impact of the long-range intra- 
band Coulomb interaction between electrons (See, e.g., 
Refs. llCHlq) . Flat kinetic energy bands, by contrast, do 
not allow screening and therefore strongly emphasize in- 
teraction effects by default. Furthermore, conventional 
perturbative treatments of the interaction (in compari- 
son to the kinetic energy) fail in flat-bands due to the 
absence of a small parameter. 

Most theoretical studies of interactions in flat bands 
use the Hubbard model with an on-site termor— The 
on-site Hubbard model incorporates just the energy 
penalty for two electrons to occupy the same site while 
ignoring the long range part of the Coulomb interac- 
tion. The on-site term leads to surprising ground states 
in the flat-band Hubbard model. For example, work by 
Nagoaka 1 - finds ferromagnetism in flat bands at specific 
fillings, near one particle per site. This is in stark con- 
trast to antiferromagnetism favored by super exchange 
in dispersive bands. 

Graphene edges, nanoribbons, and dots present phys- 
ical systems hosting flat bands. Theoretical modeling 
typically relies on the on-site Hubbard model to make 



FIG. 1: (Color online) Top: Schematic of a zig-zag nanorib- 
bon of carbon atoms. Bottom: Schematic of a ferromagnetic 
crystaf with one electron for every three unit cells in one band. 
The shaded areas correspond to single unit cells and the ar- 
rows indicate aligned electron spins. 



predictions. For example, work studying flat bands 
in on-site Hubbard models of zig-zag nanoribbon o 20 ' 21 
uses meanfield theory to argue for ferromagnetic states 
along nanoribbon edges but antiferromagnetic coupling 
between edges. An ab initio calculation^, and a work us- 
ing both the weak-coupling renormalization group and 
the density-matrix renormalization-group calculation 23 
provide similar results. 

Motivated by recent experiments on graphene 
nanoribbons^ we construct interacting lattice models 
of electrons in flat-band nanoribbons. We focus on zig- 
zag nanoribbons because here, in contrast to arm-chair 
ribbons, two flat bands arise near the Fermi level even 
in the absence of adsorbates. 3 In the top panel of Fig. [T] 
we schematically show a zig-zag nanoribbon where Rq 
(~ 2. 46 A) labels the width of a unit cell along the rib- 
bon (x direction) and L y labels the number of zig-zag 
chains across the ribbon (y direction). At low densities 
the absence of intra-band screening in flat bands sug- 
gests that the long-range part of the Coulomb interaction 
is relevant. We therefore construct models that include 
even the long-range part of the interaction. We choose 
to model a very specific regime: flat-bands in zig-zag 
nanoribbons, because we expect the absence of conven- 
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tional screening to cause flat-band electrons to order in 
a way which is completely distinct from electrons in bulk 
graphene. 

The goal of our work is to establish a set of working 
Hamiltonians of zig-zag nanoribbons. We construct a 
single-particle basis of Wannier functions. We use our 
basis to compute the interaction matrix elements. We 
then establish a projection protocol that sets up approx- 
imate flat-band models. Projection into flat bands de- 
localizes basis states due to quantum interference. The 
resulting flat-band models are highly non-trivial (incor- 
porating two bands, long-range interactions, and spin) 
and can lead to many quantum ground states even in the 
absence of significant dispersion. We make simple esti- 
mates of the low energy properties of our models at odd 
denominator fillings of a single band. 

We argue that, at low densities, the long-range part of 
the Coulomb interaction supports ferromagnetic quan- 
tum crystals (bottom panel of Fig. [IJ . Crystalline order 
projected into the flat band incorporates quantum su- 
perpositions because basis states delocalize. At low fill- 
ings direct spin exchange leads to an effective Heisenberg 
model. Our simple estimates therefore predict ferromag- 
netic crystalline order in certain parameter regimes. Our 
work sets the stage for more accurate studies of our mod- 
els with a general class of Jastrow-correlated wavefunc- 
tions that apply to flat bands <2£ 

Our protocol differs from conventional band-structure 
calculations. Flat bands, in contrast to dispersive bands, 
are, by default, strongly interacting. Conventional ap- 
plications of density functional theory accurately model 
the effect of core electrons while making very local ap- 
proximations for the Coulomb interaction between mo- 
bile electrons. Flat bands require accurate treatment of 
the long-range portion of the unscreened Coulomb inter- 
action between otherwise mobile electrons. 

In Section |TT] we consider the band structure that 
arises from non-interacting tight-binding models of 
zig-zag nanoribbons. Two flat bands are identified. In 
Section [TTT] we construct localized single-particle basis 
states, orthonormal Wannier functions, from carbon 
■n z orbitals in the honeycomb lattice model of zig-zag 
nanoribbons. Sections IIVI and [V] use the Wannier 
functions to explicitly compute Coulomb interaction 
matrix elements for one and two flat bands, respectively. 
Section IVII defines a projection scheme which limits the 
total many-body model to the flat-band portion of the 
single-particle spectrum. Section IVIII sorts terms in the 
many-body model to argue that, at low fillings, energet- 
ics favor ferromagnetic quantum crystals. Section IVIIII 
summarizes and looks forward to more accurate studies 
of the models constructed here. 



II. FLAT BANDS IN ZIG-ZAG GRAPHENE 
NANORIBBONS 

We consider interacting electrons hopping among car- 
bon sites forming zig-zag graphene nanoribbons (Fig.QJ. 
We first model the electrons in a simple non-interacting 
tight-binding picture. The single-particle tight-binding 
Hamiltonian is£ 

H o = -tJ2 (4c m + h.c), (1) 

(n,m) 

where the hopping integral is t ~ 2.7 eV for graphene 2 
and the sum is along bonds of the honeycomb lattice. 
The second-quantized operator & n creates a fermion at a 
site n. Labels n and m indicate lattice sites, in contrast 
to labels for unit cells, k, I, used in the following. 

Two bands near the Fermi level flatten for large rib- 
bon widths^ An example band structure for a narrow 
width, L y — 4, is shown in Fig. [2] Near the fermi sur- 
face, the conduction band (upper band, u) and valence 
band (lower band, d) are nearly degenerate for wavevec- 
tors q in the region qRo £ [27r/3, 47r/3] and form flat 
bands. For larger widths the bands flatten considerably. 

We examine the band width with simple ansatz flat- 
band single-particle states^ Considering states in the re- 
gion qR € [2tt/3, 47r/3] with even L y : 

<t>±(q,y) = (<f>A(q,y),±(l>B{q,y)) T 

= {{-u q )y-\±{-l)v-\u q ) L «-v) T , (2) 

for y = l,...,L y where u q = 2cos(gi?o/2), the energy 
dispersion in band T = u,d can be computed analytically: 

\E r (q)\ « \Hq,y) T Ho(<i)^y)\/\Hq,y)\ 2 

= t{l-u 2 q )u^/{l-ufy), (3) 

with 

(u q .. \ 

V .. 1 u q J 

Figure [5] compares Eq. ^ with the exact results from 
Eq. CQ). 

Eq. ([3]) can be used to determine the bandwidth. For 
partially filled lattices a narrow range of single-particle 
basis states will be occupied. The bandwidth for states 
in the flat-band sector vanishes for ribbons with large 
width: 

\E r (\q-w\^n/3)\^^, (4) 

From this estimate we see that band dispersion plays 
a small role for dilute ribbons with increasing ribbon 
widths. 

A vanishing bandwidth, due to quantum interference, 
leaves the interaction as the dominant term in the many- 
body Hamiltonian for electrons. For dilute ribbons we 
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FIG. 2: (Color online) The dot-dashed lines indicate the en- 
ergy eigenvalues of Eq. Q versus wavevector for a nanoribbon 
of width L y = 4 and a g-space mesh of N = 44. The solid 
line shows the approximate expression for the energy, Eq. ((3j ■ 
Two flat bands form near qRo — n. In the large L y limit, the 
bands flatten for 2n/3 < qRo < 4n/3. 



will work in the approximation that Hq adds an overall 
constant energy shift to the spectrum. The full Hamilto- 
nian adds the unscreened Coulomb interaction: 



total 



(5) 



In the following we treat the dispersion as a small correc- 
tion to the interacting term. We project the Hamiltonian 
into the basis of fiat-band states. Our model becomes: 



Hi 



total 



^r(<?)4xrVr 

qGBZ,<T,r 

— > constant + vt-^HyV 



FB> 



(6) 



where the first equality is written in terms of the creation 
(annihilation) operators c qo . r (c qo . r ) for Bloch states at 
wavevector q and band V in the Brillouin zone (BZ), 
which are related to the operators for single-particle basis 
states by a Fourier transform: 



— r 



(7) 



qGBZ 



Here TLj is the lattice vector of the jth unit cell, N defines 
the number of unit cells and g-space mesh, and a G {t, i} 
denotes spin. Vp B denotes projection into flat bands such 
that the many-body eigenstates are constructed from 
Bloch states with qRo € [27r/3, 47r/3]. Many-body states 
incorporating these values of q will have essentially no 
kinetic energy. We consider this model as a centerpiece 
to understanding the electronic properties of flat-band 
nanoribbons at low densities. 

To explore possible many-body states in zig-zag 
nanoribbons we construct an accurate form for Eq. ([5]) 



Upper band (u) 




Lower band (d) 




FIG. 3: (Color online) Two-dimensional Wannier functions 
plotted as a function of position in the lattice for a ribbon of 



width L-u 



4. The Wannier functions tend to localize near 



the ribbon edges. 



in the flat-band basis. We note that the absence of any 
dispersion excludes intra-band screening as in ordinary 
Fermi liquids. Thus many-body eigenstates are deter- 
mined entirely by the interplay between various terms 
in the interaction. It is therefore crucial to accurately 
determine the interacting terms in Eq. ([6]) as prescribed 
by our choice of single-particle basis. To construct an 
accurate single-particle basis we revisit the underlying 
simple tight-binding model formed from overlapping 7r z 
orbitals. We construct orthonormal Wannier functions 
from these orbitals. The Wannier functions will serve as 
single-particle basis states, allowing the construction of 
competing terms in a many-body model. 



III. SINGLE-PARTICLE BASIS STATES: 
FLAT-BAND WANNIER FUNCTIONS 

In this section we construct a set of single-particle basis 
states in nanoribbon flat bands. We superpose carbon 
7r z orbitals to form orthogonal Wannier functions. The 
Wannier functions will then, in later sections, be used to 
accurately determine interaction matrix elements. 

In an isolated band the Wannier functions are given 
by: 



W j (r) = W (r-R j ) = 



V 



d qe -^R^ q (r), (8) 



BZ 



where D is the dimension, V is the volume of unit cell. 
The Bloch functions are * q (r) = J2m=i C mqXmq(r), 
with M atomic sites per unit cell. 

To make contact with first principles calculations on 
graphene nanoribbons^ we form Bloch functions from 
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Hi. Rewriting W at the origin gives: 



27T/3 7T 47T/3 

Wo 



2n/3 n 47T/3 
W 



FIG. 4: (Color online) Left panel: Schematic of the energy 
dispersion for a wide ribbon with the Fermi level between the 
degenerate energy bands u and d. In this regime low lattice 
filling allows us to accurately ignore the finite dispersion near 
the band edges. Right panel: The same as the left panel 
but at larger fillings of the upper u band. Here the flat- 
band approximation will only be a good approximation if the 
Coulomb interaction is much larger than the band width. 



carbon tt z orbitals, 4>(r) = y^/nze ^ r . The basis states 
become Xmq {r) = (1/,/N) E^" 1 e** R '0(r - Rj - T m ), 
where T m is the location of the mth atom in the unit 
cell. 

The coefficients C mq and energy eigenvalues E(q) are 
obtained from diagonalization of the secular equation: 



0- l H{q) C q = E(q)C q , 



(9) 



where the matrix H follows from the tight-binding 
Hamiltonian H : H{q) mn = J drx^ q (r)H Xnq{r) and 
the elements of the overlap matrix O are given by 
O mn = J drxtng(r)Xnq(r)- The eigenvectors C q 
< ('■ ,,. .... ("■,:,;}' yield the coefficients used in the defi- 
nition of the Wannier functions. In the tight-binding 
approximation we set O mn proportional to the elements 
of the identity matrix, 5 rnn . 

We solve Eq. ([9]) to construct orthonormal Wannier 
functions. We consider a one-dimensional lattice of unit 
cells along the nanoribbon. The discrete wavevectors be- 
come q = (2itq/N Rq)5c. The Wannier function located 
at Ri is then: 



N-l 



(10) 



The Wannier functions defined in this way are unique for 
a, D = 1 single band model^ but for higher dimensions 
and with more bands they are not necessarily unique ^ 
We choose a specific set of single-particle basis states by 



enforcing C„ 



|Cmg| at the edge atomic site m = 1. 



As a result we obtain a set of real Wannier functions 
symmetric about the x axis. 

The above Wannier function can be written as a sum- 
mation over all local atomic orbitals (f>(r) located at sites 



M N-l 



W (r) = N f ^2 X! a ™^( r ~ r ™)' 



(11) 



with weights a m j — ^2^^ C mq e L2lTq: >/ N and normaliza- 
tion constant Nr. The coefficients a completely deter- 
mine our choice of basis. 

We can extend our calculation of the Wannier functions 
to include both the upper and lower bands. A denser 
sampling in momentum space (i.e., larger N) yields more 
accurate Wannier functions. In practice, we find that 
the Wannier function has already converged when taking 
N = 44 for L y = 4. The Wannier functions of upper 
and lower bands for the same sample ribbon are shown 
in Fig. [31 We note that the Wannier functions localize 
symmetrically about x = with an extension of less than 
four unit cells. The Wannier functions are also symmetric 
(antisymmetric) along y for the upper (lower) band. 

The flat-band Wannier functions constructed here cor- 
respond to a specific choice of single-particle basis. By 
constructing superpositions of these functions we can 
equivalently construct a model using basis states local- 
ized on either edge of the ribbon via a simple rotation 
in the two-band space. Viewed in this way our model 
implicitly includes inter-edge coupling in narrow ribbons 
because we work in the basis of u and d bands as opposed 
to a two-edge basis. 

Our approach can be used to model graphene edges. 
Our study applies to the edge states of very wide ribbons 
provided we superpose our u and d band Wannier func- 
tions to construct left and right edge Wannier functions. 
Our model can then be used to study edges of very wide 
ribbons. But we stress that our model cannot apply to 
the electrons in the center of graphene because we have 
considered bands in nanoribbons that carry over only to 
edge states in the wide ribbon limit (For a discussion 
see Ref. [3j). In what follows we focus on narrow ribbons 
and only consider Wannier functions in the u and d band 
basis. 



IV. ONE-BAND COULOMB MODEL 

Interaction effects determine the low energy properties 
of Eq. ([5]) in the absence of significant dispersion. When 
the chemical potential lies between the nearly flat bands 
of zig-zag nanoribbons, the Coulomb interaction sets the 
dominant energy scale and mitigates response. Figure 0] 
shows schematic band structures for a wide ribbon with 
the chemical potential at the band degeneracy (left) and 
far from the flat-band region (right). In what follows we 
focus on dilute systems corresponding to the left panel. 
We can, as a first approximation, assume that the valence 
band is inert and that only the conduction band, u, will 
be active under external probes. Projection into the flat 
u band implies that the Coulomb interaction alone oper- 
ates in the massively degenerate subspace formed from u 
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band single-particle basis states. In this section we will 
consider the u band only. In the following section we will 
construct a model of both the u and d bands. 

We consider an unscreened Coulomb interaction in a 
single band: 

E ^jkl4a £ ] a 'C ka ,C l(J , (12) 

where the second-quantized operators c| CT (c ia ) create 
(annihilate) a fermion with spin a in a Wannier state cen- 
tered at the ith unit cell. The matrix elements V depend 
on the basis. We can rewrite the Coulomb interaction in 
the u band in a suggestive form: 

H v = Vo E n^riji + VijUjUj - JjjSj ■ Sj 

i i<j i<j 

+ \ E Wt^A^. (is) 

Here, the single-component and total density operators 
are = c| CT c i(T and rii = + rxjj,, respectively. The 
spin operators S< = (1/2)X) CTCT ' 3L r <r . .'C. , are defined 
in terms of the Pauli matrices <r. 

Eq. (|13p keeps all terms in the full Coulomb interac- 
tion. We compute the matrix elements in the basis of 
Wannier functions in the u band. Integral equations for 
the coefficients are given in the appendix, Eqs. (|24|) . The 
first term is the ordinary single-site Hubbard term which 
is the only term that is commonly used in models of flat- 
band nanoribbons (See, e.g., Refs. I20I and[2ll). The sec- 
ond term captures the diagonal portion of the Coulomb 
interaction at long range. The absence of a dispersion 
implies that these terms can be relevant and must be 
kept in accurate models, especially at low fillings. The 
third term, the direct exchange term, favors ferromag- 
netism for > 0. The last term represents remaining 
off-diagonal terms due to the Coulomb interaction. We 
find, by direct calculation, that the last terms are very 
small compared to the other terms for a single band. 

We compute coefficients in Eq. (|T3")) explicitly. We 
perform the integrals in Eqs. (|24p by approximat- 
ing the exponential part of the n z orbital, 0(r), 
as a linear combination of three Gaussian functions: 
X; s 7 s (128/3 s 5 /7r 3 ) 1 / 4 ze~' 9sr2 . We obtain the parameters 
7 S and (3 S from the STO-3G package. 2? Data for fitting 
the 7r z orbital with £ = 1.72 are listed in Table HI For 
numerical results shown here and in the following sec- 
tions, we use the Bohr radius, ao = 0.53A, as the unit of 
length and the Coulomb energy e 2 /47reao (~ 27.2 eV in 
vacuum) as the unit of energy. 

Table |TT] lists the coefficients computed for an L y = 4 
ribbon. As we see, all coefficients are positive and can 
be sorted by Vq > > Jij > 0. The ground state can 
be determined by an interplay between leading terms in 
Eq. (|13p and the chemical potential. These coefficients 
suggest that partially filled single bands support the for- 
mation of ferromagnetic crystals. However, the large 



TABLE I: Fitting parameters for the Gaussian approximation 
to the 7r z orbital with £ = 1.72. 



s 


1 


2 


3 


Is 
Ps 


0.15591627 
2.9412494 


0.60768372 
0.6834831 


0.39195739 
0.2222899 



TABLE II: Matrix elements for one-band (u band) case for 
inter-unit cell separations of up to 4Rq. 



Fo=2.24xlO _1 




1 


2 


3 


4 


Jij 


2.34xl0~ 2 


4.68xl0 -3 


9.21xl0" 4 


1.69 xlO -4 


Vij 


1.43xl0 _1 


9.56xl0" 2 


6.83xl0" 2 


5.23xl0 -2 



Coulomb interaction may cause mixing between the u 
and d bands. In the next section we construct a two- 
band model. 



V. TWO-BAND COULOMB MODEL 

We now consider a more comprehensive two-band 
model. The u and d bands in the flat-band region 
are essentially degenerate for wide ribbon widths. The 
Coulomb interaction can in principle favor occupancy of 
both bands or the occupancy of a single band. Accurate 
estimates of coefficients in the full two-band model will 
allow exploration of the two-band energy landscape to 
determine the band occupancy in future work. 

We construct Wannier functions in both the u and d 
bands. The Hamiltonian is dominated by the following 
terms: 

Hv = XX nir t^ 




»<j,r 

+ HE ( v ij n im jr ' - J'ij s ir ■ s jr >) 

i<j r^r' 

i<j rVr' erer' 

+ v;;4 ra c] T , a ,c. T , a ,c jT(T ). (i4) 

We have checked, by direct calculation, that other terms 
involving three and four centers are much smaller than 
terms kept in Eq. (fl"4"|) . Here we see the Hubbard and 
ferromagnetic terms as in the one-band case. The last 
term indicates a non-trivial band exchange term. The 
integrals for all coefficients are listed in the Appendix. 

Eq. (UJ| presents a central result of our work. The two- 
band model must be studied for different fillings and dif- 
ferent widths to determine expected ground states. Ta- 
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TABLE III: Matrix elements for the two-band case with L y = 
4 for inter-unit cell separations of up to 4Rq. 



l/ d =2.28xl0- 1 
l4=l-91xl0 _1 








F tl =2.24xl0- 1 
J I 'i=1.32xlO~ 1 






\i - j\ 


1 
1 




o 
Z 




q 


A 
4 






v d 


1.44x10" 


-l 


Q ^1 v1H" 

y .01 x iu 


-2 


6.79xl0" 2 


5.21x10" 


-2 


l.Ul 




1.43x10" 


-l 


y .oo x iu 


-2 


6.83xl0" 2 


5.23x10" 


-2 


1 09 
1 .UZ 


v: 3 


1.46x10" 


-l 


9.56x10" 


-2 


6.81xl0~ 2 


5.22x10" 


-2 


1.02 




2.60x10" 


-2 


3.04x10" 


3 


5.75xlO" 4 


1.09x10" 


4 






2.34x10" 


-2 


4.68x10" 


-3 


9.21xl0" 4 


1.69x10" 


4 






1.62x10" 


-2 


3.14x10" 


-3 


6.54xl0 -4 


1.27x10" 


4 






2.06x10" 


-2 


7.44x10" 


3 


2.95xl0" 3 


1.35x10" 


-3 






1.05x10" 


-2 


1.82x10" 


-3 


3.49xl0" 4 


6.43x10" 


-5 





bles IIIII and IIVI show numerically computed coefficients 
for two example widths, L y — 4 and 10. 

The tables show that the electron configurations are 
determined primarily by the diagonal components of the 
Coulomb interaction (rows 1-3). These rows are nearly 
equal indicating a band symmetry, as expected. These 
rows govern the charge degrees of freedom. Rows 4-6 
govern the spin degrees of freedom. The positive elements 
support ferromagnetism. The last two rows give rise to 
band exchange effects. 

We construct a simple fitting form for the first three 
rows. We note that the coefficients V£ and V ss can be 
thought of as a softened Coulomb interaction between 
smeared charges located at separate unit cells i and j. 
For large separations the charges appear as point charges 
and interact through the Coulomb interaction but at 
short ranges our basis states smear the electron charge 
over the width of the ribbon. We approximate and 

Vjj with a convenient analytic form: 

y. . ■ « ( e2 \ a °/ R ° - (is) 

where the fitting parameter D w is dependent on the 
width of the ribbon and can be determined with a numer- 
ical fitting as shown in Figs. [5] and [6l The last column of 
Tables [TTTI and IIVI shows D w obtained by fitting. 

Eq. (TIB"]) can be used to approximate the coefficients in 
Eq. (fl4l at low filling. At low filling and determine 
the configuration of charges. It then suffices to consider 
spin exchange terms at the separations fixed by Vy and 

. We use this procedure to suggest possible low energy 
solutions to Eq. (fT4"l) . 

VI. FLAT-BAND PROJECTION 

The flat-band limit, Eq. (jB)), establishes a unique set of 
non-perturbative models. In this section we construct a 
set of operators that allow fiat-band projection of models 
constructed in the previous sections. In the following 



TABLE IV: The same as Table [TTTI but for L y = 10. 



l/ d =1.21xl0- 1 








V u =1.17xl0- 1 






v;i=i.o3xio _1 








j-^.goxio- 2 








1 




2 




3 


4 




D w 


V d 
*ij 


8.93x10" 




6.96x10" 




5.47xl0" 2 


4.45x10" 




2.09 


v. u 
"tj 


8.74x10" 




6.84x10" 




5.41xl0 -2 


4.41x10" 




2.15 


V-- 


9.08x10" 


2 


6.93x10" 


2 


5.44xl0" 2 


4.43x10" 


. 2 


2.12 


J d 

•J I] 


2.82x10" 


2 


4.80x10" 


3 


LlOxlO" 3 


3.36x10" 


-4 




T u 
J y 


2.65x10" 


2 


6.99x10" 


3 


1.51xl0" 3 


4.30x10" 


-4 




ij 


1.65x10" 


2 


4.68x10" 


3 


1.15xl0" 3 


3.10x10" 


-4 






1.75x10" 


2 


1.00x10" 


2 


5.82xl0" 3 


3.52x10" 


-3 






1.28x10" 


2 


2.85x10" 


3 


6.30xl0 -4 


1.85x10" 


-4 





section we will then use the projected models in simple 
estimates of the low energy physics. 

To enforce flat-band projection we limit all g-space 
sums to the flat-band region (FBR) qRo G [2n/3, 47r/3]. 
We can therefore project into a single band by consider- 
ing a flat-band operator that limits itself to the FBR: 

S l = ^E E e-^-^cL. (16) 

/ qeFBR 

This operator creates states centered around the unit cell 
at Rj. We note that the states created by this opera- 
tor have finite overlap with neighbors at Rj+i when the 
flat-band region does not encompass the entire Brillouin 
zone. In the limit that the flat band encompasses the en- 
tire Brillouin zone the overlap between neighboring states 
vanishes and we have St — > ff. , Thus, the projection 
into a flat band that incorporates only a fraction of the 
Brillouin zone delocalizes basis states. 

We can rewrite our model in terms of projected den- 
sity and spin operators. The single-component and to- 
tal projected density operators are pi a = b\ a b irJ and 
Pi = pi-f + pii, respectively. The projected spin oper- 
ators are defined as: 

M^E E ^'^U^, d7) 

aa > q.q'GFBR 

We stress that the projected operators do not exhibit 
ordinary commutation relations because the underlying 
operators create overlapping states, i.e., (0|6j+i&t|0) ^ 0. 

The projected Hamiltonian can be rewritten entirely 
in terms of the above projected operators. Starting from 
an unprojected model, we impose projection using the 
following replacements: c — > b,n — > p, and S — > $. For 
example, the flat-band projected Coulomb interaction in 
the u band becomes: 

KUvV,, = \,y^prp,. ■ Y^WjIKPj -S, 

i i<j i<j 

+ \ E v^A^A/k- ( 18 ) 
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FIG. 5: The diagonal component of the inter-band Coulomb 
interaction (Vjj) for a zigzag graphene nanoribbon with L y — 
4 and N = 44. Circles are from numerical evaluation of 
Eqs. dHJ. The solid line is a fit with Eq. (JTSJ and D TO = 1.02. 




FIG. 6: The same as Fig. but for L y = 10 with D w = 2.12. 

The projected two-band model can also be obtained with 
a similar replacement applied to H v d . 

VII. LOW ENERGY PROPERTIES 

We use flat-band projection to discuss possible low en- 
ergy states of Eq. ([6]) based on simple energetic argu- 
ments. A detailed quantitative analysis of low energy 
states is beyond the scope of the present work. We make 
progress by ordering terms according to dominant energy 
scales. We then focus on example lattice fillings. 

To consider low energy solutions of Eq. ([6]) we first 
examine the kinetic term. The kinetic term enforces a 
flat-band projection provided the chemical potentials lies 
near the flat band, i.e., Fig. 0^,. It is then sufficient to 



require that many-body eigenstates of Hy utilize Bloch 
states with qR G [27r/3,47r/3]. We can analyze Eq. (fl"4")) 
with this g-space restriction by using projected operators 
constructed in the previous section. 

We first point out an intrinsic energetic ordering to 
each of the terms in Eq. (fT4|) . We rewrite each of the 
terms according to an approximate ordering by energy 
and in the projected space: 

i,j',r,r' 

H - -"Band-exch j 

(19) 

where we have redefined the diagonal Coulomb terms: 

V\tf = Vl v Vl =d ' r ' =u = Vr, and v\=f = V$, oth- 
rr' 

erwise V ^ = 0. (Note that our direct calculations 

r^r' r=r' 

find Vj,- fa Vfj .) We have also redefined the off- 

t i i -rr^r' . — r=d.r'=u . 

diagonal exchange terms: J i<t ^ = J-j, J u = J^, 

and J 4< j = Jfj, otherwise J ^ = 0. The last term in 
Eq. (|T9|) corresponds to the last term in Eq. (fM)) . 

We can understand the low energy properties of the 
first three terms in Eq. (fT9]l at a few specific fillings. 
Considering an inert d band, we assume that the u band 
is partially filled at odd denominators, v u — l/(2p + 1), 
where p — 1,2,.... (y indicates the number of particles 
per basis state.) Ignoring i?Band-cxch allows a decompo- 
sition of basis states into the u and d bands. An inert d 
band implies that the inter-band interaction leads to an 
overall shift of the chemical potential. A strong external 
gate bias canceling this shift should be able to maintain 
the u-band filling v u = l/(2p+ 1). 

In the limit of commuting projected density operators 
it is well known 29 that the first terms in Eq. (TiTH) lead to a 
charge order, i.e., one-dimensional Wigner crystals with 
lattice spacing 2p + 1 . We therefore expect that the u- 
band electrons form a classical Wigner crystal in the limit 
that the flat band encompasses the entire Brillouin zone. 
The bottom panel of Fig. [1] depicts a classical crystal 
configuration in a single spin state. 

In the limit that the projected density operators do 
not commute, the case for zig-zag nanoribbons, we pre- 
dict quantum crystals in partially filled bands. Quantum 
crystals arise, in direct analogy to Wigner crystals, as 
eigenstates of the projected density operators. For ex- 
ample, a trial quantum crystal state at v u = l/(2p+ 1) 
in spin state a is given by: 

n s L- +J >«i°>' ( 2 °) 

3=0 

This trial state appears to minimize the energy of the first 
two terms in Eq. (|19p by separating flat-band charges by 
an average of 2p unit cells. Thus the first two terms 
in Eq. (|19[) impose a rigid charge order in the u band. 
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However, the charges are significantly delocalized. A fi- 
nite overlap among neighbors implies that the charges 
exist in a superposition of several different unit cells at 
once: a quantum crystal. 

Provided a rigid charge ordering we consider the 
next lowest energy scale: low energy spin properties of 
Eq. (|f 9[) . We approximate the spin-spin coupling with 
an effective Heisenberg model for the u-band particles at 
v u = l/(2p+l): 



(21) 



Eq. (l2~Tj) applies to the case of a single band at odd de- 
nominator filling. 

The ground states of Eq. (|21l) are ferromagnetic quan- 
tum crystals. The low energy spin excitations are fer- 
romagnetic magnons. The underlying rigid charge order 
enforces a large magnon wavelength. At v u = l/(2p+ 1) 
spin wave theory yields excitation energies: 



Our flat-band model, Eq. (fT9|) . sets the stage for more 
accurate analyses with a combination of numerics and 
many-body wavef unctions. The absence of a small pa- 
rameter calls for a combination of variational studies and 
diagonalization to verify proposed ground and excited 
states.— In addition to crystals discussed here, uniform 
quantum liquids are also possible^ 

The models constructed here focus on key physics of in- 
teracting flat bands but exclude several realistic effects. 
In experiments on graphene nanostructures many cor- 
rections may be required before making a detailed com- 
parison with experiment. For example, edge roughness, 
defects, and substrate disorder can destroy the flat-band 
approximation. Furthermore, inter-band screening has 
also been ignored in the current study. While intra- 
band screening was implicitly incorporated in our model, 
screening from nearby bands could lead to corrections to 
the pure Coulomb model studied here, e.g., RKKY-type 
interactions i 32 i 33 



h0J q — 2 J Qj2p . 



. 1 [l-cos((2p+l)E ?)] 



(22) 



This dispersion offers a clear indicator of ferromagnetic 
crystals in the spin degrees of freedom. 

At finite temperatures the Mermin- Wagner theorem 
asserts that spin-spin correlations decay with a finite 
length scale in the one-dimensional Heisenberg model. 30 
Thus, ferromagnetic ordering holds only up to small 
length scales. The spin-spin correlation length at non- 
zero temperatures, T, for the Heisenberg chain with ex- 
change coupling J is: 3 - 



AJ 



(2 P +1)R 4T 



1 + B(ST/J) 



,(23) 



where A w 1.1 and B « 0.65. Our results sug- 



gest that for L y 



10 at T 



IK with J? i+3 



1.5 x 10 _3 (e 2 /47re ao) ~ 473-ftT the correlation length is 
£r/(2p+ l)i?o ~ 133. Thus about 390 unit cells contain- 
ing 130 u-band electrons are included in the formation of 
a fully magnetized domain at v u = 1/3 for these param- 
eters. 



VIII. SUMMARY AND OUTLOOK 

We constructed interacting flat-band lattice models of 
zig-zag nanoribbons. A single-particle basis of orthonor- 
mal Wannier functions were built from carbon 7r z orbitals 
in a honeycomb-ribbon lattice. The single-particle basis 
was used to explicitly compute the Coulomb matrix ele- 
ments for two ribbon widths, L y — 4 and 10. The total 
model, Eqs. JS]) and (IT4")) . was then projected into the flat 
bands of the single-particle spectrum. The projected flat- 
band model, Eq. (|19[) . suggests ferromagnetic quantum 
crystal ground states. 
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X. APPENDIX 

The coefficients in Eqs. [TBI and [T4l are given by: 
d 2 rd 2 r' 



V r = 



J- - 2 
u %3 — z 



r — r 

d 2 rd 2 r' 



^|Wor(r)W r(r')r, 



r — r 



-^(r)^ r (r)^ r (r')^(r') 



^ = /^l^r(r)^r(r')| 2 -ijS, 
j;, = 2 / flf^w: u (r)W jd (r)W lu (r')W; d (r'), 



\r — r' 

^ - / 1 —x\W* l (r)W jd (r>)\>--J ij , 



|r - f , 

I T—nW: u (v)W. ld {v)W JU {r')W* d {v'), 



V- = 



|r — r i 
d 2 rd 2 r' 



-W* u (r)W JU (r)W td (r')W; d (r'), 



|r - r' 

d Yd t i i 

V,,U ^ I T v —^ W ^Wiu(v)W* u {v)W ku (r).m 



The last term is used only in Eq. (|13 
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